I'm gonna overwrite a lot of this notebook's old content. I changed the way I'm calculating wt, and wanna test that my training worked.
In [23]:
from pearce.emulator import *
from pearce.mocks import cat_dict
import numpy as np
from os import path
In [24]:
import matplotlib
#matplotlib.use('Agg')
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()
In [25]:
from GPy.models import GPKroneckerGaussianRegression
from GPy.kern import *
import h5py
In [26]:
from pearce.emulator import GPKroneckerGaussianRegressionVar
In [27]:
training_file = '/home/users/swmclau2/scratch/xi_gg_zheng07_cosmo_v3/PearceXiggCosmo.hdf5'
test_file = '/home/users/swmclau2/scratch/xi_gg_zheng07_cosmo_test_v3/PearceXiggCosmoTest.hdf5'
In [ ]:
accs = []
log_accs = []
noise = []
for rbin in xrange(18):
f = h5py.File(training_file, 'r')
Ys = []
Yerrs = []
for i in xrange(40):
Ys.append(f['cosmo_no_%02d'%i]['a_1.000']['obs'].value[:, rbin])
Yerrs.append(f['cosmo_no_%02d'%i]['a_1.000']['cov'].value[:, rbin, rbin])
n_hods = 100
start_idx = 0
X1 = f.attrs['cosmo_param_vals']
X2 = f.attrs['hod_param_vals'][start_idx:start_idx+n_hods]
Y = np.vstack(Ys)[:, start_idx:start_idx+n_hods]
Yerr = np.vstack(Yerrs)[:, start_idx:start_idx+n_hods]
noise.append([Y, Yerr])
f.close()
# how to add training errors?
K1 =RBF(input_dim=7, ARD = False)# + RBF(input_dim=7, ARD = False)#+ Linear(input_dim = 7, ARD = False) + Bias(input_dim=7)# + White(input_dim=7)
K2 = RBF(input_dim=4, ARD = True)# + RBF(input_dim=4, ARD = False)#+ Linear(input_dim = 4, ARD = False) + Bias(input_dim=4)# + White(input_dim=4)
model = GPKroneckerGaussianRegressionVar(X1, X2, Y, Yerr, K1, K2)#, noise_var = 0.01 )
model.optimize_restarts(num_restarts=10, verbose = True)
print K1.param_array
print K2.param_array
print
f2 = h5py.File(test_file, 'r')
Y2s = []
for i in xrange(35):
Y2s.append(f2['cosmo_no_%02d'%i]['a_1.000']['obs'].value[:, rbin])
testX1 = f2.attrs['cosmo_param_vals']
testX2 = f2.attrs['hod_param_vals']#[:100]
testY = np.vstack(Y2s)#[:, :100]
f2.close()
predY, _ = model.predict(testX1, testX2)
med_acc, mean_acc = np.median( np.abs( (10**predY[:,0] - 10**testY.flatten(order='F'))/10**testY.flatten(order='F') ) ), \
np.mean( np.abs( (10**predY[:,0] - 10**testY.flatten(order='F'))/10**testY.flatten(order='F') ) )
print rbin, med_acc, mean_acc
accs.append((med_acc, mean_acc))
log_accs.append((np.median( np.abs( (predY[:,0] - testY.flatten(order='F'))/testY.flatten(order='F') ) ), \
np.mean( np.abs( (predY[:,0] - testY.flatten(order='F'))/testY.flatten(order='F') ) ) ))
print
In [ ]:
accs = np.array(accs)
In [ ]:
noise = np.array(noise)
noise.shape
In [ ]:
Y, Yvar = noise[:, 0], noise[:,1]
In [ ]:
exp_Yerr = np.sqrt(Yvar)*10**(Y)
exp_accs = np.mean(exp_Yerr/(10**Y), axis = (1,2) )
In [ ]:
f = h5py.File(training_file, 'r')
#print f.attrs.keys()
r_bins = f.attrs['scale_bins']
f.close()
rpoints = (r_bins[1:]+r_bins[:-1])/2.0
In [ ]:
plt.plot(rpoints, exp_accs)
plt.loglog()
In [ ]:
plt.plot(rpoints, accs[:,0], label = 'Median Error')
plt.plot(rpoints, accs[:,1], label = 'Mean Error')
plt.plot(rpoints, exp_accs, label = 'Shot noise')
plt.ylim([0.001, 0.5])
plt.ylabel('Perc. Accuracy')
plt.xlabel(r'$r$ [Mpc]')
plt.legend(loc = 'best')
plt.loglog();
In [ ]:
log_accs = np.array(log_accs)
In [ ]:
plt.title('Log Acc')
plt.plot(rpoints, log_accs[:,0], label = 'Median Error')
plt.plot(rpoints, log_accs[:,1], label = 'Mean Error')
plt.plot(rpoints, (np.sqrt(Yvar)/np.abs(Y)).mean(axis = (1,2)), label = 'Shot noise')
plt.ylim([0.001, 0.5])
plt.ylabel('Perc. Accuracy')
plt.xlabel(r'$r$ [Mpc]')
plt.legend(loc = 'best')
plt.loglog();
In [ ]: